Neel order in square and triangular lattice Heisenberg models 
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Using examples of the square- and triangular-lattice Heisenberg models we demonstrate that the density 
matrix renormalization group method (DMRG) can be effectively used to study magnetic ordering in two- 
dimensional lattice spin models. We show that local quantities in DMRG calculations, such as the on-site 
magnetization M, should be extrapolated with the truncation error, not with its square root, as previously as- 
sumed. We also introduce convenient sequences of clusters, using cylindrical boundary conditions and pinning 
magnetic fields, which provide for rapidly converging finite-size scaling. This scaling behavior on our clusters is 
clarified using finite-size analysis of the effective a-model and finite-size spin-wave theory. The resulting greatly 
improved extrapolations allow us to determine the thermodynamic limit of M for the square lattice with an er- 
ror comparable to quantum Monte Carlo. For the triangular lattice, we verify the existence of three-sublattice 
magnetic order, and estimate the order parameter to be M — 0.205(15). 

PACS numbers: 74.45.+c,74.50.+r,71.10.Pm 



Two-dimensional (2D) quantum lattice systems studied in 
condensed matter physics can be divided into two types: those 
with a sign problem in quantum Monte Carlo (QMC), and 
those without one. This is because recent developments in 
have enabled remarkably accurate large-scale 
studies of the latter systems, such as the square-lattice Heisen- 
berg model (SLHM) [4]. In contrast, the former systems, 
such as the triangular lattice Heisenberg model (TLHM) and 
other models with geometric frustration, are often the subject 
of controversy even regarding questions of what type of or- 
der, if any, is present. For the TLHM, it is only recently 
that the rough agreement between several theoretic alf?] and 
numerical[6, i 3i methods has made a convincing case that 
the model has three-sublattice, non-collinear 120° order. 

The density matrix renormalization grouplQ] (DMRG) is 
not subject to the sign problem, it has an error which can be 
systematically decreased by keeping more states, and even 
with modest computational effort it is extremely accurate 
for one dimensional and ladder systems. For 2D systems, 
the computational effort grows exponentially with the width. 
Ameliorating this effect is the very systematic behavior of the 
DMRG results versus the number of states kept, enabling the 
use of extrapolations to improve the accuracy. The extrapo- 
lation of the energy versus the truncation error e (also known 
as the discarded weight) to the limit e ^ often can improve 
the accuracy of the energy by nearly an order of magnitude. 
For observables other than the energy, extrapolation has been 
more problematic and is much less used. 

In this Letter we show that the difficulty in extrapolating 
local measurements A is due to the incorrect assumption that 
the error A A ^ e^^^. In fact, the simplest way to measure 
local quantities within DMRG makes analytic in e. The 
resulting improved extrapolations greatly improve one's abil- 
ity to measure order parameters in two dimensional systems. 
We demonstrate this approach with a study of the SLHM and 
TLHM systems. For the SLHM, the results for the on-site 
magnetization, extrapolated in both truncation error and sys- 
tem size, are about as good as the best published OMC.fl lloll 



For the TLHM, our new results for the magnetization are com- 
parable to the best series expansionjlSi] and GFMCI^l results. 

Another limitation of DMRG is a large loss of accuracy 
if periodic boundary conditions (BCs) are used lengthwise. 
As part of our treatment, we demonstrate an approach using 
cylindrical BCs on ^ Ly clusters and pinning magnetic 
fields. We show that with an appropriate choice of the aspect 
ratio a — L^/Ly, quantities such as the staggered magneti- 
zation scale much more rapidly to the thermodynamic limit 
than in widely used methods based on correlation functions 
on Lt, = Ly clusters with periodic BCs in both directions. 



We consider S — \ Heisenberg model 
H = J ^ ^ Si ■ Sj 



(1) 



on square and triangular lattices, where (ij) denotes nearest 
neighbor sites, and we set J = 1. We consider L^. x Ly sys- 
tems with periodic BCs in the y direction, and open BCs with 
pinning in the x direction. For the SLHM we consider both 
the standard orientation of the lattice and one tilted by 45°. In 
all cases we apply a staggered pinning field corresponding to 
infinite pinning on the edges of an auxiliary {L^ + 2)xLy sys- 
tem, e.g. ±0.5 for the standard orientation SLHM. Since our 
DMRG program conserves total Sz, for the TLHM it is not 
possible to pin all three sublattices simultaneously. Instead, 
we only pin in the z direction, pinning one sublattice (point- 
ing down), with the other two free to rotate in a cone. Thus we 
expect one sublattice in large systems to exhibit (Sz) ~ ~M, 
and the other two +M /2. 

We focus on the resulting onsite magnetization Mc — 
I (Sz) I in the center column of the system. For any fixed as- 
pect ratio a = L^/Ly, Mq approaches its thermodynamic 
limit, Mq, as L^, Ly — > oo. For L^ ^ Ly, the system looks 
more one-dimensional and we expect Mq to approach Mq 
from below. For Ly ^ L^, the strong pinning dominates and 
we expect an approach from above. We utilize intermediate 
values of a to accelerate the convergence with system size. 

First we discuss the convergence of DMRG and extrapo- 
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FIG. 1: Measurements of (Sz) for a site in the middle of the cluster 
with pinning fields applied on the ends, as a function of the trunca- 
tion error e. The results are normalized by the result extrapolated to 
£ — » 0. The solid lines are quadratic fits to the data. The 6\/3 x 3 
triangular cluster, rotated 90°, is shown. The length of the arrows is 
proportional to {Sz}, and pinning fields were -0.25, -0.25, 0.5. 



lations in the truncation error e — the sum of the density ma- 
trix eigenvalues which are discarded at each step. If the trun- 
cation of density matrix states were made starting from the 
exact ground state t/jQ, then the truncation error and energy 
error would vary as (to leading order) e ~ AE ^ \Atp\'^, 
where Aijj = ip — ipo, and ip is the new approximate ground 
state 1 1 ij]. For further discussion of energy extrapolation, see 

Refs. iilillB]. For measurements of an operator A other 

than the Hamiltonian, standard variational arguments imply 
an error proportional to (At/iI^I^o)^ and thus oc e^/^. 

Consider the special situation where ip is the lowest energy 
state within an incomplete basis B. Let C be the complement 
of B. Note that ijj is an exact eigenstate in the complete basis 
of a modified Hamiltonian in which the off-diagonal terms 
connecting B and C are set to zero. Label these coupling 
terms XV, where A is an expansion parameter Assuming ^ is 
close to the true ground state, ipo, Wip is small, and one can 
consider XV as a small perturbation. The leading term in A4>, 
neglecting energy denominators, is oc XVip, which is in C. 

Now consider a change of basis for C, negating each basis 
function. This sends A —A. Since the energy is inde- 
pendent of the change of basis, E{X) is even and we expect 
analytic behavior for E{X^). For the exact ground state t/jq^ 
the change of basis switches the sign of the C coefficients. 
The truncation eiTor e is (ideally 1 11]) the sum of the squares 
of these coefficients, and is therefore also an even function of 
A. Consider an operator A which is block diagonal within the 
B /C split. Its expectation value would also be independent 
of the change of basis, and thus an analytic function of X^. 

Within DMRG, the seemingly restrictive assumption that 
the operator A is block diagonal is easily satisfied for a local 
operator, such as Sz ■ Consider one particular DMRG step, and 
consider measuring an A which acts only on one or both of the 
central two sites, not part of the truncated left and right blocks. 
As part of the DMRG step, one finds the ground state ijj within 
the current reduced basis {B). Applying A on ijj creates a 
state which is exactly represented within this basis; therefore 
A is block diagonal. At this step only a few operators can be 



measured accurately, but as the algorithm sweeps through the 
lattice all local operators can be measured. 

To utilize this analytic behavior in an extrapolation, one as- 
sumes that successive sweeps, which become increasingly ac- 
curate as the number of states kept is increased, corresponds 
to decreasing A. A better (but still approximate) description 
of the calculation is that the ground state is approached by 
taking the most significant states out of the truncated basis C 
and putting them in B, not by making A smaller We expect 
that in the limits of large numbers of states kept the two types 
of approaches are roughly equivalent. Then, both the energy 
and central-site operators should have polynomial (i.e. ana- 
lytic) dependence on the truncation error, and one can expect 
well-behaved polynomial extrapolations, i 14il 

In Fig. 1, we show the behavior of {Sz) as a function of e 
for two modest sized systems where essentially exact results 
could be obtained. The results show no signs of nonanalytic 
behavior as e ^ 0, and are fit nicely with a quadratic form. 
We have experimented to find a reliable way to extrapolate to 
e — > 0, and have adopted the following simple procedure: we 
utilize only the most accurate decade of data available, and 
fit it with a cubic polynomial. The error bars assumed for 
the purpose of the fit are proportional to e. The extrapolation 
can be checked by a fourth order fit, or a quadratic fit over 
a smaller range. If these extrapolations agree well, we take 
as a rough error estimate the empirical parameter 0.2 times 
the size of the extrapolation from the last data point. If the 
extrapolations do not agree well, we run the calculation longer 
if feasible, or raise the error estimate substantially. 

The implications of the analytic behavior in e are signifi- 
cant: local measurements for fixed e are more accurate than 
previously thought, and the extrapolation e — > improves re- 
sults substantially and provides reasonable error estimates. 

We now turn to finite size effects. Previous QMC studies of 
the magnetization M have utilized correlation functions mea- 
sured in periodic L x L systems, and extrapolation in 1/L 
for the quantity Mq. The leading term varies as 1/L with 
a substantial coefficient. The expansion in 1/L for the peri- 
odic L X L SLHM is known in detail from chiral perturbation 
theory, allowing Sandvik to determine Mq — 0.3070(3) using 
only systems up to L = 16[4J. For the TLHM, chiral pertur- 
bation results are not available, and less robust QMC methods 
must be used, making extrapolation to L cx3 much more 
difficult. For example, Capriotti et. al. extrapolated Green's 
function Monte Carlo results with M'^ > 0.13 for L < 10 down 
to M§ - 0.04 for L ^ cx) to obtain Mq = 0.205(10). Other 
estimateslSJ range as high as Afo = 0.266. 

It is known that the leading 1 /i-scaling of the order pa- 
rameter M in the 2D Heisenberg systems is universal and 
is determined by the long-wavelength spectrum of the prob- 
lem, namely by the massless spin waves. fisll We have ana- 
lyzed the effect of the aspect ratio a = L^/Ly on the scal- 
ing for pinned cylindrical and for periodic clusters using both 
finite-size scaling within an effective cr-model and the finite- 
size spin-wave theory (FSSWT). A key conclusion from both 
methods is that the coefficient in the 1 /L correction to M de- 



3 




(x-l)/(L^-l) 1/L,. {x-l)/(L^-l) 



FIG. 2: (Color online), (a) FSSWT results for the SLHM showing 
the magnetization pattern M{x) as a function of position along the 
cluster X, for two representative clusters. A/q — 0.3034 is the bulk 
value of the staggered magnetization within the SWT. (b) Mc (M) 
vs 1/Ly results for various aspect ratios a = L^/Ly by FSSWT for 
periodic BCs (upper two sets) and cylindrical BCs (lower sets). In 
Ref. [l^ M was extracted from the correlation function and differs 
from our results in the higher order (l/L^) terms. 



pends on a and, for special aspect ratios Oc, vanishes, leav- 
ing corrections of order O^i/L?). The two methods agree 
exactly on the values of Uc for nontilted and tilted square- 
lattice clusters: for periodic systems, ac = 7.0555, while 
for cylindrical systems, for Mc in the middle of the cluster, 
ttc — 1.7639, almost exactly four times smaller. The values of 
Qfc are controlled by the cluster geometry and boundary condi- 
tions through the placement of the allowed wavevectors near 
the zeros of spin-wave energy: for periodic SLHM systems, 
one has k = x^)' whereas for the cylindrical-pinned 

geometry case, k = ( j^^^^ , ) • The factor of four improve- 
ment in the aspect ratio for the latter is due to the shift by jj^^ 
away from the ordering vector. The effective-model analysis 
determines the 1/L correction term up to an unknown factor, 
but the zero crossing is independent of it. 

The FSSWT produces parameter- free, approximate results 
for M = \{Sz) \ for all sites. Fig. |2ja) shows M{x) vs x for 
two representative clusters. Due to suppression of the long- 
wavelength spin fluctuations the magnetization is enhanced 
near the boundary. The asymptotic fall-off of the magne- 
tization away from the edge can be shown to be AI{x) w 
Mo + a/x, where a^=°° = (7rV8)~^ These FSSWT results 
are in a good agreement with the DMRG data for the SLHM 
in the non-tilted clusters shown in Fig. Oa). One can see that 
already for the Lj; x 6 clusters Mc provides a good estimate of 
asymptotic 2D value Mq when the aspect ratio is near a = 2. 

Figs- 12b) and [3b) show Mc versus l/Ly for cylindrical 
BCs, obtained by the FSSWT and DMRG, respectively. Also 
shown are the results for the LxL systems with periodic BCs, 
in Fig. |2|b) by FSSWT from this work and from Ref. 1^ 
and in Fig. [3lb) by QMC using standard coiTelation function 
methods, Ref. 4. Clearly, even for the same aspect ratio, the 
finite-size effects in the cylindrical BC clusters are 3-4 times 



FIG. 3: (Color online), (a) M{x) vs x DMRG results for the 
SLHM, for different aspect ratios. The line labeled "2D" and the 
solid diamond are the QMC L ^ oo extrapolated result, Mq = 
0.3070(3)l4i]. (b) Mc vs 1/Ly results from DMRG for the SLHM. 
The upper two curves are periodic QMC a = 1 results for A/ [4] . 

smaller than in the periodic systems. The FSSWT agrees pre- 
cisely with the effective theory on the value of ac = L7639 for 
eliminating the leading 1 /L-term. This is in a good qualita- 
tive agreement with the DMRG data, but the DMRG seem to 
indicate consistently higher values of ac ~ L9. We have also 
performed QMC calculations lflTll for the SLHM with periodic 
BCs. With the largest clusters up to 20 x 160 the "magic" 
aspect ratio is ac ~ 7.5, also higher than the effective the- 
ory value 7.0555. While we cannot exclude a change in the 
behavior on larger lattice sizes, this seems to indicate some 
insufficiency of the effective theory analysis. 

In Fig. [3lb) DMRG results for Mc for lattices ranging up 
to 20 X 10 are shown. For the 20 x 10 system up to to = 2400 
states were kept, with the run taking about 40 hours single- 
core time on a 2.6 GHz Mac Pro. This yielded a truncation 
error of order 10~^, a variational energy with an estimated 
accuracy of a part in 10^, an extrapolated energy accurate to a 
few parts in 10^, and an uncertainty in Adc of about 0.0007. 

More accurate DMRG results can be obtained for 45° tilted 
latticesi 18], allowing more detailed fits. For example, on a 
32/V2 x 8V2 system, the energies and Mc were roughly 2 
times more accurate than for the 20 x 10 nontilted system, 
and the finite size effects were smaller. The improved behav- 
ior comes from how DMRG sees the width of the system (the 
number of sites on the boundary of the left or right block) ver- 
sus the physical dimension-the greater spacing by a factor of 
\/2 in the tilted case accounts for the improvement. In Fig. 
Ua) we show results for Mc versus a = L^/ Ly for various 
Ly near the value a = 1.925 where the curves nearly inter- 
sect. The intersection of such curves as — > cx) provides a 
simple determination of both ac and Mq. The resulting value 
of ac, based on the available sizes, is somewhat larger than 
that given by FSSWT and the continuum analysis. The values 
of a are discrete because we have integral lattice dimensions. 
Performing a least squares fit of this data to the expression 

Mc(a, Ly) = Mq + a(a - ac)/Ly (2) 
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FIG. 4: (Color online). DMRG results for the 45° tilted SLHM. 
(a) The solid lines are straight segments connecting the discrete data 
points from different lattice sizes, with Ly — ly\f2. The two dashed 
lines show the bounds on the QMC result. |4(] (b) A three parameter 
fit to the data from (a), as discussed in the text. 



we obtain Mq = 0.3067, = 1.9252, and a = -0.1580. In 
Fig. Ub) we show a representation of this fit. The solid lines 
are based on the fit; the data points for a = 1.9 and a = 1.925 
are obtained from linear extrapolation along the lines shown 
in (a). The result for Mq is consistent with, and of comparable 
accuracy to the best QMC result. 

For the triangular lattice, we have studied a variety of clus- 
ters and pinning fields; these results consistently supported 
that the triangular system has the three-sublattice 120° order 
found in other studies. The cluster orientation shown in Fig. 
1 seems to be the most convenient and efficent for a DMRG 
analysis to obtain AIq. Our DMRG results for comparable lat- 
tice sizes are only slightly less accurate than for the SLHM. 

Unfortunately, the finite size analysis for the TLHM is 
much less accurate. The allowed widths in the preferred ge- 
ometry must be multiples of 3, and our results for Ly — 12 
are of low accuracy, leaving only Ly — 3, 6, 9. CuiTently, we 
do not have comparable analytical guidance, such as predic- 
tions for the optimal aspect ratio, for the triangular case. In 
Fig. |5] we show results for the TLHM with this orientation 
and also for lattices rotated by 90°. The scaling behavior ap- 
pears to be quite similar to the SLHM, but with a somewhat 
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smaller ac ~ 1.6 — 1.7. Assuming this behavior, we estimate 
Mq = 0.205(15). The results for the tilted clusters seem to 
have larger finite size effects and are less useful. Our result is 
consistent with recent QMC and series expansions for Mq for 
the TLHM liil. 



In conclusion, we have developed improved techniques for 
studying ordering in 2D lattice systems using DMRG, making 
DMRG competitive with QMC and series expansion methods 
for the 2D Heisenberg model on square and triangular lat- 
tices. These include proper scaling of local quantities with 
the discarded weight, and the use of non-traditional cluster 
geometries and BCs to improve finite-size scaling. These 
latter techniques can be used with other methods besides 
DMRG. We acknowledge the support of the NSF under grant 
DMR-0605444 (SRW), and the DOE under grant DE-FG02- 
04ER46174(ALC). 



FIG. 5: (Color online). Mc versus aspect ratio for various widths for 
the TLHM, from DMRG. The two curves labeled with ly come from 
clusters rotated by 90°, with Ly = lyVS. 
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